function X = myHarris(image,thresh)
% will return a vector X containing the coords of the corners
sigma = 1;%blurring factor
% coefficients for derivation, given in:   Hany Farid and Eero Simoncelli 
%"Differentiation of Discrete Multi-Dimensional Signals" IEEE Trans. Image 
%Processing. 13(4): 496-508 (2004)
p = [0.037659  0.249153  0.426375  0.249153  0.037659];
d1 =[0.109604  0.276691  0.000000 -0.276691 -0.109604];
Iy = conv2(d1, p, image, 'same');
Ix = conv2(p, d1, image, 'same');%actual differentiation

%gaussian blur:
h = fspecial('gaussian',[3 3],sigma);
Iy2 = imfilter(Iy.^2,h);
Ix2 = imfilter(Ix.^2,h);
Ixy = imfilter(Ix.*Iy,h);

%improved harris corner measure, corner strength image!
Hm = (Ix2.*Iy2 - Ixy.^2)./(Ix2 + Iy2 + eps);%eps is a small numerical error 

% Extract local maxima by performing a grey scale morphological
% dilation and then finding points in the corner strength image that
% match the dilated image and are also greater than the threshold.
radius = 1;
sze = 2*radius+1;                   % Size of dilation mask.
Hm_dl = ordfilt2(Hm,sze^2,ones(sze)); % Grey-scale dilate.

% Make mask to exclude points within radius of the image boundary.
bordermask = zeros(size(Hm));
bordermask(radius+1:end-radius, radius+1:end-radius) = 1;

% Find maxima, threshold, and apply bordermask
X_im = (Hm==Hm_dl) & (Hm>thresh) & bordermask;
[c,r] = find(X_im);                % get row,col coords.

X = [r c];

